home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 4: GNU Archives / Linux Cubed Series 4 - GNU Archives.iso / gnu / glibc-1.09 / glibc-1 / glibc-1.09.1 / sysdeps / generic / pow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-31  |  6.6 KB  |  216 lines
  1. /*
  2.  * Copyright (c) 1985, 1993
  3.  *    The Regents of the University of California.  All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)pow.c    8.1 (Berkeley) 6/4/93";
  36. #endif /* not lint */
  37.  
  38. /* POW(X,Y)  
  39.  * RETURN X**Y 
  40.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  41.  * CODED IN C BY K.C. NG, 1/8/85; 
  42.  * REVISED BY K.C. NG on 7/10/85.
  43.  * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
  44.  * Required system supported functions:
  45.  *      scalb(x,n)      
  46.  *      logb(x)         
  47.  *    copysign(x,y)    
  48.  *    finite(x)    
  49.  *    drem(x,y)
  50.  *
  51.  * Required kernel functions:
  52.  *    exp__D(a,c)            exp(a + c) for |a| << |c|
  53.  *    struct d_double dlog(x)        r.a + r.b, |r.b| < |r.a|
  54.  *
  55.  * Method
  56.  *    1. Compute and return log(x) in three pieces:
  57.  *        log(x) = n*ln2 + hi + lo,
  58.  *       where n is an integer.
  59.  *    2. Perform y*log(x) by simulating muti-precision arithmetic and 
  60.  *       return the answer in three pieces:
  61.  *        y*log(x) = m*ln2 + hi + lo,
  62.  *       where m is an integer.
  63.  *    3. Return x**y = exp(y*log(x))
  64.  *        = 2^m * ( exp(hi+lo) ).
  65.  *
  66.  * Special cases:
  67.  *    (anything) ** 0  is 1 ;
  68.  *    (anything) ** 1  is itself;
  69.  *    (anything) ** NaN is NaN;
  70.  *    NaN ** (anything except 0) is NaN;
  71.  *    +(anything > 1) ** +INF is +INF;
  72.  *    -(anything > 1) ** +INF is NaN;
  73.  *    +-(anything > 1) ** -INF is +0;
  74.  *    +-(anything < 1) ** +INF is +0;
  75.  *    +(anything < 1) ** -INF is +INF;
  76.  *    -(anything < 1) ** -INF is NaN;
  77.  *    +-1 ** +-INF is NaN and signal INVALID;
  78.  *    +0 ** +(anything except 0, NaN)  is +0;
  79.  *    -0 ** +(anything except 0, NaN, odd integer)  is +0;
  80.  *    +0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
  81.  *    -0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
  82.  *    -0 ** (odd integer) = -( +0 ** (odd integer) );
  83.  *    +INF ** +(anything except 0,NaN) is +INF;
  84.  *    +INF ** -(anything except 0,NaN) is +0;
  85.  *    -INF ** (odd integer) = -( +INF ** (odd integer) );
  86.  *    -INF ** (even integer) = ( +INF ** (even integer) );
  87.  *    -INF ** -(anything except integer,NaN) is NaN with signal;
  88.  *    -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
  89.  *    -(anything except 0) ** (non-integer) is NaN with signal;
  90.  *
  91.  * Accuracy:
  92.  *    pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
  93.  *    and a Zilog Z8000,
  94.  *            pow(integer,integer)
  95.  *    always returns the correct integer provided it is representable.
  96.  *    In a test run with 100,000 random arguments with 0 < x, y < 20.0
  97.  *    on a VAX, the maximum observed error was 1.79 ulps (units in the 
  98.  *    last place).
  99.  *
  100.  * Constants :
  101.  * The hexadecimal values are the intended ones for the following constants.
  102.  * The decimal values may be used, provided that the compiler will convert
  103.  * from decimal to binary accurately enough to produce the hexadecimal values
  104.  * shown.
  105.  */
  106.  
  107. #include <errno.h>
  108. #include <math.h>
  109.  
  110. #include "mathimpl.h"
  111.  
  112. #if (defined(vax) || defined(tahoe))
  113. #define TRUNC(x)    x = (double) (float) x
  114. #define _IEEE        0
  115. #else
  116. #define _IEEE        1
  117. #define endian        (((*(int *) &one)) ? 1 : 0)
  118. #define TRUNC(x)     *(((int *) &x)+endian) &= 0xf8000000
  119. #define infnan(x)    0.0
  120. #endif        /* vax or tahoe */
  121.  
  122. const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;
  123.  
  124. static double pow_P __P((double, double));
  125.  
  126. double pow(x,y)      
  127. double x,y;
  128. {
  129.     double t;
  130.     if (y==zero)
  131.         return (one);
  132.     else if (y==one || (_IEEE && x != x))
  133.         return (x);        /* if x is NaN or y=1 */
  134.     else if (_IEEE && y!=y)        /* if y is NaN */
  135.         return (y);
  136.     else if (!finite(y))        /* if y is INF */
  137.         if ((t=fabs(x))==one)    /* +-1 ** +-INF is NaN */
  138.             return (y - y);
  139.         else if (t>one)
  140.             return ((y<0)? zero : ((x<zero)? y-y : y));
  141.         else
  142.             return ((y>0)? zero : ((x<0)? y-y : -y));
  143.     else if (y==two)
  144.         return (x*x);
  145.     else if (y==negone)
  146.         return (one/x);
  147.     /* x > 0, x == +0 */
  148.     else if (copysign(one, x) == one)
  149.         return (pow_P(x, y));
  150.  
  151.     /* sign(x)= -1 */
  152.     /* if y is an even integer */
  153.     else if ( (t=drem(y,two)) == zero)
  154.         return (pow_P(-x, y));
  155.  
  156.     /* if y is an odd integer */
  157.     else if (copysign(t,one) == one)
  158.         return (-pow_P(-x, y));
  159.  
  160.     /* Henceforth y is not an integer */
  161.     else if (x==zero)    /* x is -0 */
  162.         return ((y>zero)? -x : one/(-x));
  163.     else if (_IEEE)
  164.         return (zero/zero);
  165.     else
  166.         return (infnan(EDOM));
  167. }
  168. /* kernel function for x >= 0 */
  169. static double
  170. #ifdef _ANSI_SOURCE
  171. pow_P(double x, double y)
  172. #else
  173. pow_P(x, y) double x, y;
  174. #endif
  175. {
  176.     struct Double s, t, __log__D();
  177.     double  __exp__D(), huge = 1e300, tiny = 1e-300;
  178.  
  179.     if (x == zero)
  180.         if (y > zero)
  181.             return (zero);
  182.         else if (_IEEE)
  183.             return (huge*huge);
  184.         else
  185.             return (infnan(ERANGE));
  186.     if (x == one)
  187.         return (one);
  188.     if (!finite(x))
  189.         if (y < zero)
  190.             return (zero);
  191.         else if (_IEEE)
  192.             return (huge*huge);
  193.         else
  194.             return (infnan(ERANGE));
  195.     if (y >= 7e18)        /* infinity */
  196.         if (x < 1)
  197.             return(tiny*tiny);
  198.         else if (_IEEE)
  199.             return (huge*huge);
  200.         else
  201.             return (infnan(ERANGE));
  202.  
  203.     /* Return exp(y*log(x)), using simulated extended */
  204.     /* precision for the log and the multiply.      */
  205.  
  206.     s = __log__D(x);
  207.     t.a = y;
  208.     TRUNC(t.a);
  209.     t.b = y - t.a;
  210.     t.b = s.b*y + t.b*s.a;
  211.     t.a *= s.a;
  212.     s.a = t.a + t.b;
  213.     s.b = (t.a - s.a) + t.b;
  214.     return (__exp__D(s.a, s.b));
  215. }
  216.